library(tidyverse)
library(janitor)
library(sf)
library(tmap)
… medmindre man ikke aner, hvordan man gør!
sfDen typiske måde at lege med spatialt data på er vha. GIS (Geographic Information Systems)-værktøjer designet til det
Programmer som R er dog løbende blevet udvidet med pakker, der gør det muligt at klare alting i det samme stykke software, som man bruger til andre ting
Det er dobbelt smart, fordi man har alting ét sted og bygget op omkring kode, der kan ændres og opdateres
Med sf er det blevet smooth sailing. Den måde, pakken håndterer det spatiale aspekt af et datasæt gør, at det ligner alle andre datasæt til forveksling
I kan kende næsten alle sf-funktioner på deres st_*-prefix (fx st_cast())
sfData ifølge {sp}
Spatialt data opbevares typisk enten som
x, y-koordinatergeojsonshapefileI kan læse geojson-filer med `sf::read_sf``:
url <-
"https://api.dataforsyningen.dk/kommuner?format=geojson"
kommuner_raw <-
read_sf(url)
rm(kommuner_raw)
Til shapefiles foretrækker (plejer?) jeg at bruge sf::st_read(). Her angiver dsn den mappe, jeres (flere!) filer ligger i. layer angiver navnet, de forskellige filer har til fælles (dvs. alt undtagen endelser, fx .shp)
europe <- st_read(dsn = "data/europe",
layer = "europe")
## Reading layer `europe' from data source
## `C:\Users\jevi\OneDrive - Epinion\Documents\GitHub\crashcourse_geoviz\data\europe'
## using driver `ESRI Shapefile'
## Simple feature collection with 48 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -10.39023 ymin: 36.02593 xmax: 41 ymax: 71.17769
## Geodetic CRS: WGS 84
I kan gemme jeres shapefiles ved hjælp af sf::st_write():
st_write(europe,
dsn = "data/europe",
layer = "europe",
encoding = "UTF-8",
delete_layer = TRUE,
factorsAsCharacter = TRUE,
driver = "ESRI Shapefile")
## Deleting layer `europe' using driver `ESRI Shapefile'
## Writing layer `europe' to data source `data/europe' using driver `ESRI Shapefile'
## Writing 48 features with 3 fields and geometry type Multi Polygon.
df$geometryDet smarte ved sf’s måde at opbevare spatialt data på er, at ALT det spatiale ligger i den liste, der hedder geometry. Den er også “fredet”, og man kan ikke slette den med select(). Hvis man vil af med den (og det vil man af og til, det vender vi tilbage til!), er der en specifik måde til det.
Udover geometry ligner spatiale datasæt en tibble som I kender dem. Bemærk, at sf tilføjer geografisk metadata, når I printer datasættet (inkl. CRS, datatype og geografisk omfang, bounding box)
europe
## Simple feature collection with 48 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -10.39023 ymin: 36.02593 xmax: 41 ymax: 71.17769
## Geodetic CRS: WGS 84
## First 10 features:
## cntry unit sub_re geometry
## 1 Vatican Vatican Southern Europe MULTIPOLYGON (((12.43838 41...
## 2 United Kingdom Jersey Northern Europe MULTIPOLYGON (((-2.082227 4...
## 3 United Kingdom Guernsey Northern Europe MULTIPOLYGON (((-2.520898 4...
## 4 United Kingdom Isle of Man Northern Europe MULTIPOLYGON (((-4.392285 5...
## 5 United Kingdom United Kingdom Northern Europe MULTIPOLYGON (((-2.539355 5...
## 6 Ukraine Ukraine Eastern Europe MULTIPOLYGON (((38.20586 47...
## 7 Switzerland Switzerland Western Europe MULTIPOLYGON (((9.35 47.598...
## 8 Sweden Sweden Northern Europe MULTIPOLYGON (((18.95645 57...
## 9 Spain Spain Southern Europe MULTIPOLYGON (((1.592676 38...
## 10 Slovakia Slovakia Eastern Europe MULTIPOLYGON (((22.47305 49...
tribble(~name,~n_champs,~lat, ~lon,
"Vejle Boldklub",5,55.71372762470475, 9.556222451831141,
"Aarhus Gymnastikforening",5,56.13193948661876, 10.196519237914524,
"AC Horsens",0,55.87162048058711, 9.857707600274574)
## # A tibble: 3 x 4
## name n_champs lat lon
## <chr> <dbl> <dbl> <dbl>
## 1 Vejle Boldklub 5 55.7 9.56
## 2 Aarhus Gymnastikforening 5 56.1 10.2
## 3 AC Horsens 0 55.9 9.86
df <- tribble(~name,~n_champs,~lat, ~lon,
"Vejle Boldklub",5,55.71372762470475, 9.556222451831141,
"Aarhus Gymnastikforening",5,56.13193948661876, 10.196519237914524,
"AC Horsens",0,55.87162048058711, 9.857707600274574)
df <- df %>%
st_as_sf(coords = c("lon", "lat"),
crs = 4326)
tmap_mode("view")
tm_shape(df) +
tm_bubbles(size = 2,
col = "n_champs")
pointsurbana <- st_read(dsn = "data/urbana",
layer = "urbana")
class(urbana$geometry)
## [1] "sfc_POINT" "sfc"
ggplot() +
geom_sf(data = urbana,
size = .01) +
theme_void()
linesrivers <- st_read(dsn = "data/rivers",
layer = "rivers")
class(rivers$geometry)
## [1] "sfc_MULTILINESTRING" "sfc"
ggplot() +
geom_sf(data = rivers) +
theme_void()
polygonseurope <- st_read(dsn = "data/europe",
layer = "europe")
class(europe$geometry)
## [1] "sfc_MULTIPOLYGON" "sfc"
ggplot() +
geom_sf(data = europe) +
theme_void()
ggplot2Vi vil gerne vise størrelsen af urbane centre i Europa
Vi har også en idé om, at mange større byer ligger tæt på floder, fordi det historisk set har været en fordel at være tæt på vand, fordi det gjorde handel lettere.
Lad os starte med at vise kortet over Europa!
Fordi vi bruger sf-pakken til det hele, kan vi bruge tidyverse-funktioner til at manipulere data med – og ggplot2 til at visualisere det!
Vi har ikke brug for at vise koordinater hen ad akserne, så vi bruger theme_void():
ggplot() +
geom_sf(data = europe) +
theme_void()
Ikke dumt! Lad os tilføje vores data over byer. Det gør vi let ved at tilføje endnu en “geometry” til vores `ggplot:
ggplot() +
geom_sf(data = europe) +
geom_sf(data = urbana) +
theme_void()
Aha! Det er jo meget fedt, at vi har data på byer rundt i hele verden, men lige her gør det nok mere skade end gavn.
Husk, at tidyverse-funktioner kan bruges på sf-datasæt, akkurat som på klassiske datasæt. Hvis der var en regionsvariabel i bydatasættet, kunne vi derfor let filtrere alt uden for Europa væk vha. dplyr::filter(). Det er der desværre ikke, så vi må ty til andre tricks.
Heldigvis har sf-pakken en lang række værktøjer, der hjælper os med at manipulere data. Fx kan vi bruge st_intersection() til at beskære ét datasæt med et andet:
eu_cities <- urbana %>%
st_intersection(.,
europe)
ggplot() +
geom_sf(data = europe) +
geom_sf(data = eu_cities,
size = .8) +
theme_void()
Meget bedre! Men, det bliver en smule tætpakket med alle de prikker. Lad os lege lidt med det visuelle:
ggplot() +
geom_sf(data = europe, fill = "white") +
geom_sf(data = eu_cities,
size = 1,
alpha = .5) +
theme_void()
Da vi beskar bydatasættet med shapefilen med Europa-kortet, tilføjede den automatisk data fra sidstnævnte.
eu_cities
## Simple feature collection with 3424 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -9.251242 ymin: 36.1835 xmax: 40.98619 ymax: 69.74173
## Geodetic CRS: WGS 84
## First 10 features:
## scalerank featurecla area_sqkm min_zoom tmp cntry unit
## 3294 9 Urban area 8.936 7.7 3294 United Kingdom Jersey
## 3293 8 Urban area 34.616 7.6 3293 United Kingdom Guernsey
## 3231 7 Urban area 18.935 7.0 3231 United Kingdom Isle of Man
## 3198 9 Urban area 10.900 7.7 3198 United Kingdom United Kingdom
## 3199 9 Urban area 11.659 7.7 3199 United Kingdom United Kingdom
## 3201 9 Urban area 10.329 7.7 3201 United Kingdom United Kingdom
## 3202 7 Urban area 80.259 7.0 3202 United Kingdom United Kingdom
## 3203 9 Urban area 6.911 7.7 3203 United Kingdom United Kingdom
## 3204 8 Urban area 32.645 7.6 3204 United Kingdom United Kingdom
## 3205 9 Urban area 10.634 7.7 3205 United Kingdom United Kingdom
## sub_re geometry
## 3294 Northern Europe POINT (-2.098328 49.19537)
## 3293 Northern Europe POINT (-2.550666 49.47685)
## 3231 Northern Europe POINT (-4.48212 54.16456)
## 3198 Northern Europe POINT (-3.312512 57.65061)
## 3199 Northern Europe POINT (-3.635281 57.62275)
## 3201 Northern Europe POINT (-4.213206 57.47909)
## 3202 Northern Europe POINT (-2.128716 57.1617)
## 3203 Northern Europe POINT (-2.467238 56.73593)
## 3204 Northern Europe POINT (-2.607162 56.59364)
## 3205 Northern Europe POINT (-2.738971 56.52185)
Det kan vi bruge til at tilføje flere farver til vores kort!
ggplot() +
geom_sf(data = europe, fill = "white") +
geom_sf(data = eu_cities,
aes(fill = sub_re),
shape = 21,
size = 1,
alpha = .5) +
theme_void()
Ikke dumt! Men det er stadig svært at udlede så meget substantielt. Lad os prøve noget nyt: bruge de data, vi fik smidt på byerne før til at summere nogle mere aggregerede mål. På den måde kan vi visualisere landene (polygoner) med information fra byerne (punkter).
Bydatasættet har ikke data på population, men der er et mål for areal i kvadratmeter. Lad os antage, at det er en fin proxy for befolkning (her antager vi konstant befolkningstæthed, hvilket selvfølgelig ikke er korrekt, men lad os køre med det).
Her har vi to spatiale datasæt:
Det kan selvfølgelig kombineres! Igen, vi kan bruge alle de tidyverse-tricks vi kender og elsker.
Fordi vi i sidste ende gerne vil visualisere polygonerne, har vi ikke brug for det spatiale element i bydatasættet. Det skal bare merges igen. Derfor sletter vi geometry-listen. Hvorfor ?
sf vil ikke merge to spatiale datasæt med en klassisk *_join(). Den vil meget hellere joine på det spatiale, men det har vi ikke brug for hergeometry’en, når I ikke skal bruge den. Lad mig vise det nedenfor:library(tictoc)
# aggreger vores areal-variabel inden for lande med spatialt data
tic()
city_cntry <- eu_cities %>%
group_by(cntry) %>%
summarize(area = sum(area_sqkm)) %>%
ungroup()
toc()
## 1.8 sec elapsed
# aggreger vores areal-variabel inden for lande med ikke-spatialt data
temp <- eu_cities %>%
st_drop_geometry()
tic()
city_cntry <- temp %>%
group_by(cntry) %>%
summarize(area = sum(area_sqkm)) %>%
ungroup()
toc() ; rm(temp)
## 0 sec elapsed
Forskellen er jo ikke enorm, men den er klar! Hvis I arbejder med større datasæt, kan det hurtigt blive en enorm fordel at tænke over, hvornår det spatiale aspekt af data mest bare er en hæmsko.
Nå, I digress. Tilbage til visualiseringen. Det nye city_cntry data kan vi let joine med vores landepolygoner med en simpel left_join(), og så er vi back in business!
df_europe <- europe %>%
group_by(cntry) %>%
summarize()
df_europe <- df_europe %>%
tidylog::left_join(.,
city_cntry,
by = "cntry")
Det gik jo smooth! Vi kan dog se, at tre lande (i polygondatasættet) ikke blev merged med noget data. Det er sådan noget man altid lige skal tjekke (og derfor jeg ELSKER {tidylog}-pakken), men her er det bare fordi der ikke var nogle større byer. Lad os for god ordens skyld bruge den information til at erstatte NA med 0 i vores nye variabel:
df_europe <- df_europe %>%
mutate(area = ifelse(is.na(area),
0,
area))
Nu kan vi lave vores nye kort!
ggplot() +
geom_sf(data = df_europe,
aes(fill = area)) +
theme_void()
Bum! Det var lige det, vi gerne ville have. Men faktisk er det lidt dumt, at Rusland er med på kortet, når det ikke rigtigt er en del af Europa – og fordi vi alligevel kun viser en lille del af det. Lad os fjerne det med en simpel filter():
df_europe <- df_europe %>%
filter(cntry != "Russia")
ggplot() +
geom_sf(data = df_europe,
aes(fill = area)) +
theme_void()
Sådan! Nu vi er i gang, kan vi godt prøve at gøre det lidt pænt.
ggplot() +
geom_sf(data = df_europe,
aes(fill = area),
color = "white",
size = .01) +
scale_fill_viridis_c(name = "Sum of Area (km^2)\n ofUrban Areas") +
theme_void() +
guides(fill = guide_colorbar(direction = "vertical",
barheight = 12,
barwidth = .5))
Det eneste, der måske kan være lidt misvisende er, at vores befolkningsmål ikke er korrigeret for landenes areal. Det er heldigvis nemt at implementere med st_area() fra sf-pakken og lidt datarwrangling:
df_europe <- df_europe %>%
mutate(poly_area = st_area(.)) %>%
mutate(poly_area = unclass(poly_area))
df_europe <- df_europe %>%
mutate(relative_area = area/poly_area)
ggplot() +
geom_sf(data = df_europe,
aes(fill = relative_area),
color = "white",
size = .01) +
scale_fill_viridis_c(name = "Urban Area as Share\n of Total Area") +
theme_void() +
guides(fill = guide_colorbar(direction = "vertical",
barheight = 12,
barwidth = .5))
Det tegner jo et lidt andet mønster.
Lad os, for at vende tilbage til den oprindelige analyse, se om datasættet over floder kan bidrage lidt!
# beskær data til kun at dække Europa
rivers <- rivers %>%
st_intersection(.,
df_europe)
ggplot() +
geom_sf(data = df_europe,
aes(fill = relative_area),
color = "white",
size = .01) +
scale_fill_viridis_c(name = "Urban Area as Share\n of Total Area") +
geom_sf(data = rivers,
color = "orange") +
theme_void() +
guides(fill = guide_colorbar(direction = "vertical",
barheight = 12,
barwidth = .5))
ggsave(plot = last_plot(),
filename = "plots/europe_static.png")
tmapTil de fleste formål vil ggplot2 være bedre til at visualisere jeres spatiale data på. Syntaxen er også noget mindre logisk, selvom den overordnede struktur minder om ggplot2’s.
tmap har dog en central fordel: det er eminent til at lave interaktive kort!
Af samme grund bruger jeg det flittigt, når jeg arbejder med/manipulerer data, fordi det gør det nemt at inspicere datakilder og -kvalitet:
# sæt til interaktiv visning
tmap_mode("view")
tm_shape(df_europe) +
tm_polygons(col = "relative_area", alpha = 0.85,
border.col = "white", lwd = .8) +
tm_shape(rivers) +
tm_lines(col = "steelblue")
Bemærk i øvrigt, at i modsætning til geom_sf() i ggplot2, skal tmap have specificeret typen af spatial data (med andre ord, hvilken geometry, der er tale om).
html <- tm_shape(df_europe) +
tm_polygons(col = "relative_area", alpha = 0.85,
border.col = "white", lwd = .8) +
tm_shape(rivers) +
tm_lines(col = "steelblue")
tmap::tmap_save(html,
"plots/europe_interactive.html")
“All models are wrong, but some are useful”
– George Box
Vi vil ofte gerne vise en globe i et 2D-rum. Problemet er, at er fysisk umuligt. Det kan simpelthen ikke lade sig gøre. Det betyder, at alle kort I nogensinde har set (glober undtaget) er forkerte på én eller anden måde. Tror I mig ikke? Se her!.
Det er her, CRS’er (Coordinate Reference Systems) kommer ind i billedet. For at “løse” det helt fundamentale problem (og fordi glober er ret besværlige at gå rundt med) er der udviklet en lang række projektioner. Vi bruger med andre ord disse matematiske projektioner til at projicere en globe (oftest vores egen klode) i et 2D-rum på en måde, der fungerer til formålet.
Alle disse projektioner er “forkerte” på en eller anden måde, og der er altid en eller anden forvrængning af størrelse, form eller retninger. Men, de er jo ikke udviklet for sjov. De er udviklet, fordi vi har brug for dem, og fordi de er brugbare. De er ikke allesammen lige brugbare, og de er sjældent lige brugbare til samme formål. Med en let omskrivning af ovenstående citat kan vi beskrive dem som:
→ All projections are wrong, but some are useful
Det nok mest berømte eksempel på en projektion (og på, at de aldrig er perfekte…) er Mercator-projektionen (EPSG 4326 / 3857, se mere her). Den er berømt og berygtet for at repræsentere visse dele af kloden (især Afrika og Sydamerika) som værende mindre, end de faktisk er. Det er jo uden tvivl rigtigt nok. Se billedet nedenfor, hvor cirklerne angiver areal relativt til faktisk areal. Den mere ukvalificerede del af kritikken går så langt som til at fremhæve, at den er “forkert”. Det er den, men det er alle alternativer også.
Der er også spekuleret meget i, at projektionen skulle være udviklet med henblik på mentalt at (re)producere europæisk dominans, fordi den fremstiller Afrika som langt mindre relativt til Europa end det faktisk er:
“These maps of Africa, drawn up by a small group of western cartographers, symbolically reinforced Europeans’ sense of control over their mapped territories and subjects, but they didn’t betray much in the way of real information. Though they would have been seen as objective and impartial at the time, in retrospect it is clear how subjective, ideologically driven, and, in many ways, fantastical they were.”
– James Wan
Det er noget sludder. Det er også et ulogisk argument. I det 15. århundrede, hvor Mercator-projektionen blev udviklet, havde europæerne ikke behov for at fremstille det afrikanske argument som værende “småt” for at legitimere kolonialisering. Der var masser af “legitimitet” at hente qua religion, “civilisation”, vanvittige idéer om “raceforskelle” osv. Tværtimod havde europæiske magthavere, om noget, interesse i at fremstille deres oversøiske besiddelser som værende så store og imponerende som muligt.
Mercator-projektionen er uden tvivl relateret til europæisk kolonialisme, men på en lidt mere indirekte måde. Den var nemlig helt eminent at navigere efter, især på meget lange maritime strækninger og med simpel teknologi. Som, fx, på turene over Atlanterhavet i det 15. århundrede. På disse længere ture giver Mercator-projektionen ikke helt den korteste rute, men det giver den mest bullet-proof kurs.
Det vigtigste er, at I bruger den samme CRS i alle jeres objekter. Et lille trick er at definere en værdi i starten, som I løbende kan referere til. På den måde skal I kun ændre den ét sted, hvis det bliver aktuelt
my_crs <- 4326
# df <- read_sf() %>%
# st_transform(4326)
#
# map <- read_sf() %>%
# st_transform(3359)
#
# st_crs(df) == st_crs(map)
#
# df <- df %>%
# st_transform(crs = my_crs)
#
# map <- map %>%
# st_transform(crs = my_crs)
#
# st_crs(df) == st_crs(map)
For mere information om brug af CRS’er i R, se:
rgdal::make_EPSG()